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ABSTRACT 

This paper considers the impact of Lyman-a coupling and X-ray heating on the 21-cm 
brightness-temperature one-point statistics (as predicted by semi-numerical simulations). The 
X-ray production efficiency is varied over four orders of magnitude and the hardness of the 
X-ray spectrum is varied from that predicted for high-mass X-ray binaries, to the softer spec¬ 
trum expected from the hot inter-stellar medium. We find peaks in the redshift evolution of 
both the variance and skewness associated with the efficiency of X-ray production. The am¬ 
plitude of the variance is also sensitive to the hardness of the X-ray SED. We find that the 
relative timing of the coupling and heating phases can be inferred from the redshift extent of 
a plateau that connects a peak in the variance’s evolution associated with Lyman-a coupling 
to the heating peak. Importantly, we find that late X-ray heating would seriously hamper our 
ability to constrain reionization with the variance. Late X-ray heating also qualitatively alters 
the evolution of the skewness, providing a clean way to constrain such models. If foregrounds 
can be removed, we find that LOFAR, MWA and PAPER could constrain reionization and 
late X-ray heating models with the variance. We find that HERA and SKA (phase 1) will be 
able to constrain both reionization and heating by measuring the variance using foreground- 
avoidance techniques. If foregrounds can be removed they will also be able to constrain the 
nature of Lyman-a coupling. 

Key words: dark ages, reionization, first stars - intergalactic medium - methods: statistical 
- cosmology: theory. 


1 INTRODUCTION 

The cosmic dark ages, during which the only source of radia¬ 
tion was the adiabatically cooling cosmic microwave background 
(CMB), ended when the first stars formed (see |Loeb & Furlan-[ 
e tto|2013) for an overview of this process). The exact nature of the 
first stars and galaxies is uncertain, but the radiation they emitted 
will have dramatically altered the state of the intergalactic medium 
(IGM). Neutral hydrogen (Hi) dominates the composition of the 
IGM until the epoch of reionization (EoR); as such it is hoped that 
the impact of these galaxies will be detectable with the 21-cm line, 
a spectral line produced by a hyperfine transition in HI ( |Field|1958| 
|1959j >. The spin temperature (T s ), defines the distribution of the 
electron population over the singlet and triplet hyperfine levels in¬ 
volved in the 21-cm transition i |Ewen & Purcel 1|[T 951 f If this is 
in equilibrium with the temperature of the CMB (T CM b) then the 
21-cm signal will not be detectable. However, radiation from stars 
breaks this equilibrium, leading to an observable signal in absorp¬ 
tion where T s < Tcmb, and in emission when T s > T CMB . 


* Email: c. watkinsonl 1@imperial .ac.uk 


1.1 The impact of the first stars 

As well as affecting the spin temperature, radiation from the first 
stars began ionizing neutral hydrogen. Most important to this dis¬ 
cussion is the production of Lyman-a, X-ray and ultra-violet radi¬ 
ation. 

(i) Wouthuysen-Field (Lyman-a) coupling:. The populations 
of the 21 -cm hyperfine levels are mixed by repeated absorption and 
re-emission of Lyman-a radiationJ^This couples the spin tempera¬ 
ture to the Lyman-a colour temperature T a . Repeated scattering of 
Lyman-a photons off atoms couples T a to the kinetic gas tempera¬ 
ture (T k ) and so T s ~ T k i)Wouthuysen| 1 952||Field| 1 958||Pritchard| 
|& Furl anetto 2006). WF coupling therefore produces fluctuations in 
the 21-cm signal, the observation of which would provide insight 
into the nature of Lyman-a sources. 

(ii) X-ray heating: An abundance of X-rays are produced by 
accretion onto compact objects, such as black holes and neutron 
stars, as well as by hot gas in the interstellar medium. These X- 
rays induce photo-ionizations resulting in primary and secondary 
electrons. It is unlikely that X-ray photo-ionizations are efficient 

1 Lyman-n photons contribute to the Lyman-a radiation field through cas¬ 
cades. 


© 2015 RAS 
















2 C. A. Watkinson & J. R. Pritchard 


enough to be solely responsible for the reionization of the universe 
( |Dijkstra et al.|2004||Hickox & Markevitch|2007{|McQuinn|2012fr ; 

however, once the IGM has been ionized to a few percent, the 
photo-ejected electrons deposit the majority (~ 65%) of their en- 


ergy as heat in the IGM ( Shull 1979; 

Shull & van Steenberg|1985, 

Furlanetto & Johnson Stoever/2010 

). Because the Wouthuysen- 


Field (WF) coupling described in (i) is likely to have started at 
an early stage, the onset of X-ray production will raise the 21-cm 
spin temperature. Observation of 21-cm fluctuations produced by 
the heating process will therefore provide insight into the nature of 
X-ray sources. 

(iii) Reionization: As the Universe evolves, Ultra-Violet (UV) 
radiation from a growing number of galaxies begins to ionize the 
IGM. These ionizing photons have a short mean-free path and so 
carve out well defined ionized regions around galaxies in an other¬ 
wise mostly neutral IGM. These grow with time until they eventu¬ 
ally merge and the Universe is reionized. The associated depletion 
of neutral gas (and thus 21-cm signal) produces 21-cm fluctuations 
whose observation will provide insight into the process of reion¬ 
ization (i.e. the nature of the sources responsible, and of the IGM) 
Whilst sub-dominant to UV ionizations. X-ray induced ionizations 
will also impact on the ionization state of the IGM (for example, 
|Qh|200l||Venkatesan et al.|2001| and |Mesinger et al.|2013) . 

The exact timing of processes (i)-(iii), and the degree to which 
they overlap, are uncertain, especially for the coupling and heating 
processes. This uncertainty depends on the nature of the stars that 
drive WF coupling and the efficiency (and relative timing) of X- 
ray production (Mesinger et al. 201 3[ >. If the X-ray efficiency is 
low enough, 21-cm fluctuations induced by X-ray heating may well 
persist into the EoR. 

Despite such uncertainties, it is expected that processes (i) - 
(iii) occurred in the order described. As the main source of X-rays is 
thought to be accretion onto compact objects, the production of X- 
rays is likely to be delayed by a few million years relative to the first 
stars igniting. Lyman-a production is coincident with formation of 
the first stars, and the emissivity to achieve Lyman-a coupling is 
much lower than that required of X-rays to substantially heat the 
IGM. Therefore Lyman-a coupling will at least commence prior to 
the onset of X-ray heating i jChen et~ah 2004). Because the mean 
free path of UV photons is very small, UV-driven reionization will 
inevitably be delayed relative to WF coupling and X-ray heating. 


1.2 Constraints on X-ray production from post-EoR 
redshifts 

The only constraints we have on the nature of high-redshift X-ray 
production come from lower-redshift observations. The dominant 
X-ray sources observed are active galactic nuclei (AGN) and high- 
mass X-ray binaries (HMXBs). X-ray emission from the hot inter¬ 
stellar medium (ISM) is also found to contribute significantly to the 
soft X-ray emission of nearby galaxies (e.g. jMineo et al.|2012b) . 

Observations of the unresolved cosmic X-ray background 
point towards AGN being the dominant contributor in the local 
universe ( [Moretti et al.||20l2) . However, the AGN number den¬ 
sity rapidly decreases at z > 3 (see |Fan et ah|[200l| and refer¬ 
ences therein), although some scope remains in the faint end of 
the luminosity function for low-mass ‘mini-quasars’ to contribute 
at higher redshifts (for example, Madau et al. 2004; Volo nteri &| 
|Gnedin|2009) . 

HMXBs are expected to be dominant at high redshift because: 
(1) in the absence of an AGN they dominate X-ray output in low- 


redshift galaxies l |Fabbiano|2006| l, and (2) their abundance is pro¬ 
portional to star-formation rate and ‘starburst’ galaxies are seen to 
increase with redshift (for example, Gilfanov et al. 2004, Mirabel 
|et al.|20nj |Mineo et al.|2012a) . Theoretical modelling also sug¬ 
gests that a high fraction of the first stars formed in binaries or mul¬ 
tiple systems would evolve into HMXBs ( |Turk et al.|2009[|Stacy| 
|et al.|2010| >. 


The X-ray spectral energy distribution (SED) and the associ¬ 
ated mean free path of X-rays determine the scale of 21-cm fluc¬ 
tuations produced by inhomogeneous heating. The SED can be fit 
with a power law where the specific luminosity Lx is proportional 
to the frequency v as Lx oc v~ a with spectral energy inde^of 
a. A harder X-ray SED, with a < 1, produces higher energy X- 
rays, which have a longer mean free path. Shocks in supernovae and 
AGN have an energy index a ~ 1.7 (Tozzi et al. 2006; McQuinn 
2012), the hot ISM produces an SED described by a ~ 3 iPacucci 


et al. 


2014), and HMXBs have a hard spectra described by a ~ 0.7 


1 (Rephaeli et al.|1995[|Swartz et al.|2004[|Mineo et al.|2012a) . 

Even if we knew the appropriate spectral energy index to describe 
the SED at high redshift, there is the matter of the luminosity’s nor¬ 
malisation. We can normalise the luminosity at high redshifts to 
match the low redshift observations, however it could be orders of 
magnitudes apart from what we observe today. 


1.3 Observing and understanding the 21-cm signal 


Given that the X-ray properties of high redshift sources are so un¬ 
certain there is much to be gained from observing these epochs. The 
first generation of 21-cm radio telescope such as lofaf(5 mwaQ 
and PAPEF0 aim to constrain reionization statistically and are al¬ 
ready starting to set interesting upper limits (see Paciga et al. 201 1[ 
|Dillon et al.|2013||Ali et al.|2015| and |Pober et al.|2015| >. However, 
next-generation instruments HER/^j and S K will not only dra¬ 
matically improve constraints on reionization, but also aim to probe 
the pre-reionization era. Telescopes seeking to measure the global 
average of the 21-cm signal, such as EDGE^Jand DARI0 should 
also provide valuable (and complementary) constraints on the EoR 
and pre-reionization epochs (Bowm an & Rogers|2010| Burns et al. 

|20l2| ). 

Observing the 21-cm line will clearly be rewarding, however 
it will be challenging as the signal is small (~ 10 mK) and fore- 


grounds will be orders of magnitude larger ( Shaver et al. 1999 

Dt 

Matteo et al. 2002b, Oh & Mack 2003| Di Matteo et al. 2004a 

). It 


is hoped that by exploiting the spectral smoothness of foregrounds 
they may be removed (e.g. Wang et al. 2006, Liu & Tegmarkj201Tl 


Paciga et al.|201 1| Petrovic & Qh|201 l||Chapman et al.|2012| Cho| 


et al.|2012[ Shaw et al.|2014| >. Alternatively, we could avoid fore¬ 

grounds by exploiting the existence of a wedge feature in the k±-k y 


2 This is related to the photon index as T = a + 1. 

3 The LOw Frequency ARray http : / /www. lof ar . org/ 

4 The Murchison Wide-held Array http://www.mwatelescope . 
org/ 

° The Precision Array to Probe Epoch of Reionization http: //eor. 
berkeley.edu/ 

u The Hydrogen Epoch of Reionization Array http:// 
reionization.org/ 

' The Square Kilometre Array http: / /www. skat ele scope .org/ 

8 The Experiment to Detect the Global EoR Signal http://www. 
haystack.mit.edu/ast/arrays/Edges/ 

a The Dark Ages Radio Explorer http://lunar.colorado.edu/ 
dare/ 
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cylindrically-binned 2D power spectrum to which foregrounds are 


confined (e.g. Datta et al. 2010 Vedantham et al. 

2012; Morales 

et al. 2012, Thyagarajan et al. 2013; Hazelton et al. 

2013 ; Liu et al. 


2014b i W It is not yet clear how well foregrounds can be miti 


gated (see |Liu et afj|2QI4a| , so it is vital that we have a strong 
understanding of the statistics of 21-cm fluctuations, even in light 
of next-generation instruments. There are also a large number of 
astrophysical parameters (many of which are degenerate with each 
other, see for example |Greig & Mesinger||2015| and |Pober et al.| 
|2015j > for which we have no constraints on in the high-redshift Uni¬ 
verse. Therefore we must also fully investigate all possibilities for 
the range of physics we might observe. 


The paper is structured as follows: In Section [2] we describe 
our simulations; in Section [TT] we study the evolution of the vari¬ 
ance and skewness during the epochs of WF coupling and X-ray 
heating; we then look at the impact of X-ray processes on the mo¬ 
ments during reionization in Section [3~2| in Section [4] we consider 
the observational prospects for constraining the nature of WF cou¬ 
pling, X-ray heating, and reionization using the moments; finally in 
Sections[5]and[6]we discuss caveats of our work and conclude. 


2 SIMULATION OVERVIEW 


1.4 Overview of this work 


Much attention has been focussed on measuring the 21-cm power 
spectrum, which has been shown to be rich with information (for 


example,!Furlanetto et al. 2004a, Zaldarriaga et al. 2004, Mellema 

et al. 

2006, Lidz et al. 2008, Pritchard & Loeb 2008; Santos et al. 

2008 

Mesinger et al.|2011| Friedrich et al. 2011; Mesinger et al. 

2013 

Sobacchi & Mesinger[2014, Greig & Mesinger|20151. How- 


ever, given the challenging nature of these observations it is also 
worth considering one-point statistics. One-point statistics have 
been shown to be information rich, are simpler to interpret, and 


will be differently affected by foregrounds (e.g. Furlanetto et al. 

2004b Wyithe & Morales 2007 Harker et al. 2009 

Ichikawa et al. 

2010, Watkinson & Pritchard 2014, Watkinson et a 

. 2015). 


The sensitivity of the 21-cm one-point statistics to coupling 
and heating has not yet been studied in detail. In this paper we in¬ 
vestigate the sensitivity of these statistics to the X-ray efficiency 
and spectral index using semi-numerical simulations. In doing so, 
we lift the assumption that T s T CM b (which is often made 
when simulating reionization) to study the impact of different X- 
ray properties on the 21-cm moments during the EoR. 

We note that during the writing of this paper a similar work 
by |Shimabukuro et aI7| ( |2014) > was submitted to MNRAS. Our work 
differs in that our simulated boxes are bigger (theirs are 200 Mpc, 
ours are 600 Mpc). This is of particular importance in studying X- 
rays because of their long mean free path, which can be up to hun¬ 
dreds of Mpc (see |McQuinn||20 i~2) . We also include peculiar ve¬ 
locities in our simulations. Our paper includes several additional 
elements: 

• We present detailed analysis of the impact of X-ray processes 
on reionization. In particular, we show that the evolution of the 
skewness is altered in the case of late X-ray heating, providing a 
useful feature for constraining such models. 

• The impact of different values for the spectral index is stud¬ 
ied, identifying a degeneracy between X-ray efficiency and spec¬ 
tral index (because both alter the amplitude of the 21-cm variance). 
This degeneracy may be broken through observations of the 21-cm 
skewness. This provides sensitivity to the X-ray spectral hardness. 

• Finally we consider the prospects for constraining X-ray 
source properties with current and future generations of radio tele¬ 
scope. In particular, we establish that even if foreground removal is 
not possible, the 21-cm variance can be accurately measured using 
foreground avoidance techniques. 


We use the latest public release version of 21CMFAST (v 1.12) for 
this work. For details of this simulation we refer the reader to 
|Mesinger & Furla netto (20071 and Mesinger et al. (20111; however 
for convenience we will summarise the main points. 

The code uses the Zel’dovich Approximation ( |ZeTdovich| 
|1970| applied to linear-density fields to generate evolved density 
(5) and velocity (dry/dr) fields. The excursion-set formalism of 
|Furlanetto et al.| ( [2004b| > can then be applied to the evolved den¬ 
sity fields to generate neutral-fraction (rc H i) fields to model UV- 
driven ionizations. The offset of the brightness temperatur^] rel¬ 
ative to that of the CMB (6Tb) can then be calculated (assuming 
T s » Tcmb) according to, 


6T h = T \ TcMB (l- e -^o), 


: 27 


1+z 
T s — Tc MB 
T s 


£Chi (1 + 6 ) 


H(z)/(l + z) 


dry / dr 


/ 1 + z 0.15 
\ 10 D m /i 2 


1/2 


D b /r\ 

0.023 ) 


ml(. 


( 1 ) 


The cosmological parameters in this equation are the Hubble pa¬ 
rameter H(z), the matter D m , and baryon density parameters fib 
(where Hi = pi /p c and p c is the critical density required for flat 
universe). 

Throughout this paper we adopt a standard ACDM cosmol¬ 
ogy as constrained by Planck, i.e. a g = 0.829, h = 0.673, 
D m = 0.315, D a = 0.685, D b = 0.049 and n s = 0.96 ( [Planck] 
[Collaboration et al.|2014| >. Our simulations produce 3D boxes with 
side L = 600 Mpc and cubic pixels of side ( p j x = 1 Mpc (initial 
conditions are calculated at double this resolution). All lengths are 
co-moving unless otherwise stated. 


2.1 X-ray heating and ionizations 

This paper focuses on the effects of X-ray heating, as such it is 
essential that the details of spin temperature fluctuations are in¬ 
cluded. 21CMFAST calculates the spin temperature according to 
[Fieldl < fl958l > ; |Hirata| < [2006] l , 


10 The wedge results from the spectral smoothness of the foregrounds 
(which one might only expect to observe at low kj_) combined with the 
chromatic nature of a radio telescope (i.e. at different frequencies the in¬ 
strument probes different scales) which smear foregrounds into larger fcy. 


11 In radio observations, the radiation intensity /„ is described in terms 
of brightness temperature, T b , defined such that /„ = B(T b ); B(T) is 
the Planck black-body spectrum and is well approximated by the Rayleigh- 
Jeans formula at the frequencies relevant to reionization studies. 
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T’cmb + x a'T a 1 + x c T k 1 

1 Xa. Xc 

Jcmb + x aT a 1 

1 X OL 


( 2 ) 


The WF-coupling coefficient is defined as x a = 1.7 x 10 11 (1 + 
z)~ l S a J a pcm _2 s _2 Hz^ 1 sr _1 , where S a ~ 1 is a factor cor¬ 
recting for detailed atomic physics and T a is the Lyman-a colour 
temperature. The second equality of Equation[2]follows because the 
collisional coupling coefficient x c « 0 during the epoch relevant 
to this work. Both S a and T a are calculated according to |Hirata| 
( |2006) l. For the purpose of our analysis, we define an ‘effective’ 
coupling co-efficient as x e ff = x a /(t + x a) assuming S a = 1. 

The Lyman-a flux ( J a ) has two main origins, (1) X-ray excita¬ 
tion (via photo-electrons) of neutral hydrogen and (2) direct stellar 
emission between the Lyman limit and the Lyman-a line. Lyman- 
n photons are redshifted into Lyman-a resonance; therefore, J a is 
calculated by integrating the photon contribution (as produced via 
mechanisms 1 and 2) over a series of concentric spherical redshift 
shells surrounding each pixel. For more detail than provided here, 
we refer the interested reader to the 21CMFAST literature listed at 
the beginning of this section as well as Flirata (2006); Pritchard & 
|Furlanetto|(2006 ; and |Pritchard & Furlanetto| ( |2007| >. 

Where WF coupling is saturated, x a T CM b and T a - « 
T~ x \ in this regime, the kinetic temperature will be tightly coupled 
to T a due to the repeated scattering of the Lyman-a photons by 
hydrogen atoms < jField| 1958) . The kinetic temperature (outside of 
FI II regions at position x and at redshift z) is calculated by tracking 
the heating history for that position. This can be calculated using 
the evolution of T^(x, z') which is predicted by, 


dr k (a:, z') 
d z' 


dt 


3fc B (l + x e ) dz 

Tk dx e 
1 + x e dz' 




Cp + 


221 drib 
3rib dz' 


(3) 


is required. The X-ray luminosity of sources is assumed to be well 
described by a power law, i.e. //emitted oc ( v/vf)~ a (where a is 
the spectral index discussed in Section [772] ) and hvo is the lowest 
energy X-ray that can escape into the IGM). 21CMFAST assumes 
that the total X-ray emission rate per second (dN/ dz”) from a 
spherical shell bounded by the redshift interval z” to z" + dz" is 
the product of the number of X-ray photons per solar mass in stars 
£x (the X-ray efficiency parameter) and the star formation rate in 
that shell (SFR,L) (i.e. dN/ dz" = (xSFRz//)0The arrival rate 
at position x and redshift z' from sources between z" and z" + dz" 
is then, 

dcj>{x, is, z', z") __ d N -if is \ -“- 1 / 1 + z" \ -“- 1 _ Tx 
dz" d z" aU ° \is 0 ) Vl + z'J 

(5) 


where r x is the optical depth of X-rays. In calculating the mean free 
path of X-rays, fluctuations of the IGM state are ignored. Note that 
this is very inaccurate during the advanced stages of reionization, 
when there are large ionized regions in an otherwise neutral IGM. 
With this approximation, the heating rate due to X-rays at position 
x and redshift z' is calculated as, 


r oo 

e x (x,z') = d v 

J Max[i/o,i/ T -i] 


X ^2(hu - F^/heat/iZiCTi 
i 


X 



dz" 


d </>/ dz" 

A-KTp 


(6) 


where r p is the proper (null geodesic) separation of z' and z". Un¬ 
der the same assumption, the ionization rate due to X-rays may be 
described by, 


(x,z) = 


j 

J1V 


Max[i/o,F T - 




d (j)/ dz" 
4nrp ’ 


(7) 


where 


Fi =(his - Ef 1 ) 


/ion,Hi /ion,HeI /ion,Hell 

inth + rsth + c;th 
-^Hl ^Hel & HeII 


+ 1 • ( 8 ) 


dx e (x,z') _ dt r 
dz' ~dz' L i( 


aACxinbfii] ■ 


(4) 


In calculating the kinetic temperature, 21CMFAST must also cal¬ 
culate the local ionized fraction x e {x,z'), which depends on the 
total baryon number density rib, the ionization rate per baryon 
Fion, the case-A recombination co-efficient oa, the clumping fac¬ 
tor C =< n 2 > / < n H > 2 (where n H describes the hydrogen 
number density and we set C = 2), and finally the hydrogen num¬ 
ber fraction / H . In addition, the kinetic temperature depends on the 
heating rate (e p ) per baryorj^jfor process p (for our discussion the 
dominant process is X-ray heating) and the Boltzmann constant k B . 

To calculate X-ray heating and ionization rates at redshift z' 
one must integrate over the full range of frequencies for which pho¬ 
tons can contribute energy to these processes. Furthermore, to ac¬ 
count for redshifted photons, another integral over redshift (z" 0 


12 The heating rate (e p ) has units erg s 1 . 

13 Implicit is the assumption that z" > z'. 


In Equations [6] to [8] a sum is taken over the species i =H I, 
He I, He II; Xi is the ionization fraction for the specie^] fi is the 
species number fraction, ct; is the ionization cross section, and 
is the species ionization threshold energy. The fraction of the pri¬ 
mary electron’s energy that is transferred to heat and secondary 
ionizations is described by /heat,i and fion,i respectively for each 
species. The unity term in Fi accounts for primary ionizations of 
species i. The heating and ionization rate are simplified further (to 
speed up the calculation) by assuming that no photons with an op¬ 
tical depth r < 1 are absorbed by the IGM and all photons with 
r > 1 are jMesinger & Furlanetto|2007]|Mesinger et al.|2011[>. 


14 The star formation rate SFR ; ,// = Pcrito^b/*[l + 

(z")] AV/ dz" df con (z",R")/ dt, where dV(z") is the co¬ 
moving volume element at z ", (z") is the PT-evolved density field 

smoothed on scale R ", /* is the fraction of baryons converted to stars, 
/coll is the collapsed fraction (which is calculated as described in Section 
12.21, and p cr it,o is the critical density at z = 0. 

13 Xi = 1 — x e for H I and He I; Xi = x e for He II. 
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2.2 UV Ionizations 

Ionizations by UV photons are calculated independently from, and 
following the simulation of X-ray heating and WF coupling. The 
code smooths iteratively around each pixel in the box from a maxi¬ 
mum radius 7? max [^]down to pixel scales _R p i x . At each smoothing 
step, on scale R, the condition / co n(x, z, R) > f© 1 is evaluated, 
if met the central pixel is marked as ionized; if this condition is 
never met then the pixel is partially ionized, accounting for both 
UV ionizations, calculated as = Cuv/coii(r’, z, R p i x ), and par¬ 
tial ionizations due to X-rays (calculated using Equation[4]l. 

The collapsed fraction / co n is calculated using the predic¬ 
tion of the extended Press-Schechter formalism ( |Bond et al.|1991| 
|Lacey & Cole|l993] l, with the minimum mass corresponding to a 
virial temperature of 10 4 K (necessary for cooling by atomic hy¬ 
drogen to be effective). The collapse fraction is normalised so that 
its mean agrees with that predicted by the parametrically fit mass 
function of jjenkins et al.j j2001 f . 

The ionizing efficiency of stars is defined as £ uv = 
/esc/*(V 7 / 6 (l+n r ec) _1 , where / eac is the fraction of ionizing pho¬ 
tons that escape, iV 7 /i, denotes the number of ionizing photons pro¬ 
duced per baryon in stars, and finally n lec is the typical number of 
times that hydrogen will have recombined. The fraction of baryons 
converted to stars /* also impacts upon the estimation of SFR used 
in the X-ray heating and ionization rates. We set £ uv = 16 for con¬ 
sistency with our previous publications, and so that the 50% point 
of reionization falls in the redshift range to which first generation 
instruments are most sensitive, whilst agreeing with observational 
constraints on the EoR. Our reionization model is thus optimistic. 

Because UV photons have a short mean free path, it is as¬ 
sumed that they will carve out large ionized regions in a mostly 
neutral IGM; it is useful then to consider the volume filling factor 
of these ionized Hll regions Q H ii. The average ionized fraction of 
the box, taking into account the X-ray ionizations discussed above, 
is Xion — Qhu + (1 — Quu^Xe. 


3 RESULTS 

In this work, we study the properties of X-ray sources by vary¬ 
ing the luminosity normalisation and spectral index. The nor¬ 
malisation of the luminosity is parametrised through the X-ray 
efficiency parameter £ x in 21CMFAST . We simulate £ x = 
[10 55 ,10 56 ,10 57 ,10 58 ], with Cx = 10 57 (roughly 1.7 X-ray pho¬ 
tons per stellar baryon) as our fiducial model. This choice is con¬ 
sistent with low-redshift constraints on the total X-ray luminosity 
per unit of star formation (with /* = 0.1). We then consider the 
hardness of the X-ray background by varying the spectral index, 
assuming values ranging from a = 0.8 (to approximate the spec¬ 
trum produced by HMXB ) to a = 3.0 (typical of a soft X-ray back¬ 
ground as produced by the hot ISM). We set a = 1.5 in our fiducial 
model, as it is in the middle of the plausible range of values for this 
parameter and is representative of X-ray hardness in the local Uni¬ 
verse (see the discussion in Section O- Unless otherwise stated 
results are from maps that have been smoothed and re-sampled to 
produce pixels with side 4 Mpc. This is to overcome the impact of 
a discretisation effect (that occurs through the creation of the non- 


16 We set Umax = 30 Mpc based on the ionizing UV photon mean free 
path at the redshifts of interest, see|Storrie-Lombardi et al. (1994); Miralda- 
|Escude||2003};|Choudhury et al.||2008jl 




Brightness temperature (mK) 


Figure 1. Redshift evolution of the spin (T s ), kinetic (Tk), and CMB (T CM b ) 
temperatures (top), and brightness-temperature PDF (bottom) at five key 
phases of its evolution for the fiducial model with f x = 10 5 ' and a = 1.5. 
All PDFs are from maps smoothed to 4 Mpc. Arrows on the top plot corre¬ 
spond to the redshifts for which we present PDFs. The nature and timing of 
these features are model dependent; in some extreme models it is possible 
that not all of the five phases we describe will exist, or that their PDFs will 
look quite different. 


linear density fields) on the moments (see [Watkinson & Pritchard] 
|2014| for details). 

Ignoring fluctuations in peculiar velocities and at a fixed red¬ 
shift (and cosmology), the drivers of brightness-temperature fluc¬ 
tuations are the density, neutral-fraction (iEhi = 1 — Xion) and spin- 
temperature fields; specifically, 

5T h oc /tZhiA , 

with (i = l — T CM b/T s and A = 1 + S. At early times, before 
the epoch of reionization, x HI ~ 1 and fluctuations in fi and A 
dominate the signal. Therefore any evolution of the brightness tem¬ 
perature away from that of the density field will be due to correla¬ 
tions between (i and A. As such we can gain insight by looking at 
the cross averages of these quantities, which can be broken up as 
follows, 

{(J, *HlA} = XHI - HI^ - + (&Shi) 

or when ami = 1: 



( 10 ) 
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z=30.62, (Xi on ) = 0 z-25.52, (x ]0 „) = 0 z=22.24, (x ion ) = 0 z=18.49, (x ion )=0.01 z=16.08, (x to „)=0.02 



Figure 2. Brightness-temperature maps of our fiducial (f x = 10 57 , a = 1.5) simulation. Maps are from boxes of side L = 600 Mpc smoothed to have 
pixels of (4 Mpc) 3 and presented with the pixel depth flattened into the page. Maps correspond to the five key phases we describe in the text; i.e. z ~ 30 
where T s ~ T c mb (left), z ~ 25 where T s < T c mb (middle-left), z ~ 22 where T s is at its minimum (middle), z ~ 18 where T s ~ T c mb (middle-right) 
and z ~ 16 where T s 3> T c mb (right)- The amplitude of the signal is model dependent, the more efficient the X-ray production the less the amplitude in 
absorption. If the X-ray efficiency is high enough, it is even possible that the signal will never enter a phase of absorption as seen in this fiducial model. 


3.1 Coupling and heating epochs 


Before considering the evolution of the 21-cm moments, we can 
build some insight by looking at the probability density function 
(PDF) of the brightness temperature. The top plot of Fig. [T] shows 
the redshift evolution of the three temperature components rele¬ 
vant to (STi,: the kinetic temperature (T^\ blue-dashed line), the 
CMB temperature (T CM b; black-dashed line w/triangles) and the 
spin temperature (T a ; red solid line). The bottom plot of the same 
figure shows the shape of the 5Tb PDF at five important phases of 
the brightness-temperature’s evolution for our fiducial model (see 
Shimabukuro et al. 2014 for discussion of the 1 — T C mb/T b PDF, 


:: mb / ’( 

)0 


which agrees with the interpretation we present below)[ 

For reference, the associated brightness-temperature maps are 
presented in Fig.[2]along with the redshift evolution of spin temper¬ 
ature, ionized fraction and x e a in Fig. [3] (top, middle and bottom 
respectively). 


• T s ~ Tcmb> z ~ 30 (Blue dot-dashed line in Fig.[l]bottom): 
WF coupling begins almost immediately in our simulations and is 
positively correlated with the density field (this can be seen in the 
plot of the cross average of the ‘effective’ WF coupling coefficient 
with density as a function of redshift in the bottom of Fig. [4}. As 
such the spin temperature in overdense regions (near the sources 
of Lyman-a radiation) is becoming coupled to the gas temperature 
(which is cooling adiabatically with the expansion of the universe). 
As a result the mean spin temperature drops below that of the CMB. 
This process produces a negatively skewed brightness-temperature 
PDF which is quite sharply peaked with the weight of the distribu¬ 
tion towards more negative brightness temperatures. At this point, 
brightness-temperature fluctuations are dominated by fluctuations 
in the density and Lyman-a flux. 

• T s < Tcmb, z ~ 25 (black dashed line with triangles in 
Fig.[T]bottom): The Lyman-a coupling coefficient and the density 
field are most strongly correlated around this epoch for all models 
presented in this paper (again refer to Fig.[4]bottom). As the Lyman- 
a coupling becomes more effective the spin temperature starts to 
evolve more rapidly towards gas temperature, and the skewness of 
the PDF becomes less negative (as the statistics of the density field 


1 ' We choose to plot the PDFs with a log y-axis as we find it better for 
visualizing skewness in the distributions. 
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Figure 3. Redshift evolution of: spin (thick lines) and kinetic (thin lines) 
temperatures (top); ionized fraction (middle); and effective coupling coef¬ 
ficient (bottom), a = 1.5 for all models. Only in the extremely X-ray- 
efficient Cx = 10 5s model do X-ray ionizations contribute substantially 
to reionization. All other models exhibit a very similar, UV-driven, reion¬ 
ization history. However, the evolution of the spin temperature is different 
between all models. 
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become increasingly influential). The variance is increasing during 
this phase. 

• T s at its minimum, z ~ 22 (red solid line in Fig.^bottom): 
Eventually the spin temperature reaches a minimum just before 
coupling fully with the (now increasing) gas temperature. From 
the PDF we can see that despite the average brightness tempera¬ 
ture being at its minimum, some more extreme pixels are already 
in emission; i.e. coupling and X-ray heating are both very strong 
in some pixels. At this point, the PDF has a positive skewness, pri¬ 
marily driven by fluctuations in the X-ray heating but amplified by 
fluctuations in the WF-coupling. This is because a region that is 
less strongly coupled will have a spin temperature closer to that of 
the CMB; a region that is both strongly coupled and more heated 
than the mean will also result in a spin temperature closer to that of 
the CMB. 

• T s ~ T C mb again, z ~ 18 (Green dashed line w/circles in 
Fig. |T| bottom): The spin temperature is now fully coupled to the 
gas temperature and is thus increasing due to X-ray heating. At this 
point, fluctuations in the X-ray heating are dominating those of the 
brightness temperature. The average brightness temperature is zero, 
and fluctuations produce a relatively even distribution of pixels in 
emission and absorption; therefore the skewness is close to zero. 
The variance is also decreasing as X-ray heating is becoming more 
homogeneous. 

• T s T cmb, 2 ~ 16 (Pink dotted line with stars in 

Fig. □ bottom): Eventually the spin temperature becomes much 
greater than the CMB temperature and heating fluctuations become 
unimportant. This results in a nearly Gaussian distribution as the 
brightness-temperature fluctuations are governed nearly entirely by 
those of the density field. Reionization by UV photons is just be¬ 
coming effective around this time. An earlier reionization model 
and/or less efficient X-ray production could mean that this Gaus¬ 
sian phase never occurs; instead there may be a phase in which 
fluctuations in both the heating and ionization fields occur at the 
same time (as seen in the extreme ( x = 10 55 , which we describe at 
length in Section [T2| i. 

It is important to note that the PDFs described are from our 
fiducial (Cx = 10 57 , a = 1.5) model. Thus, these five points may 
be observed at different redshifts; the evolution of the PDFs will 
also vary quantitatively in different models. Furthermore, if X-ray 
production is either extremely efficient, or extremely inefficient, 
then the evolution of the various temperatures and therefore the 
PDFs will be qualitatively different from the fiducial model. 

3.1.1 Efficiency of X-ray production 

Fig. a (top) shows the redshift evolution of the brightness- 
temperature PDF’s variance. The variance is zero at very high 
redshift for all models. It then increases with decreasing redshift, 
driven by a slight positive correlation between the density field and 
Tf ; i.e. the spin temperature is smaller in overdense regions, be¬ 
cause WF-coupling is strongest in the vicinity of sources and dur¬ 
ing this phase Tk < T c mb. This is illustrated by the evolution of 
{TcmbS/T s ) shown in the middle plot of Fig. [7] The evolution of 
the variance plateaus briefly as the average spin temperature drops 
towards the average gas temperature (although note this is less ev¬ 
ident in the £ x = 10 58 as X-ray heating occurs so early). Eventu¬ 
ally an anti-correlation between the density field and T a _1 develops. 
By this point, WF-coupling fluctuations are minimal (see the bot¬ 
tom plot of Fig. [4]l and so this effect is caused by the underdense 
regions being less heated by X-rays than those closer to sources; 



Ffed shift 


Figure 4. Redshift evolution of: brightness-temperature variance (top) 
which exhibits two peaks prior to that of reionization; the first is driven by 
(bottom) and the second by ((Tqmb /T s )S) (middle). Thin lines on 
the middle plot correspond to the evolution of the ionized fraction for each 
model. Statistics are calculated from maps resolved to 4Mpc, and a = 1.5 
for all models. 

i.e. the spin temperature is smallest in underdense regions where 
there are less X-ray sources. In all but the (j x = 10 55 model, the 
variance is largest when this anti-correlation is maximized. As we 
will see, the £ x = 10 55 model enters this phase during the early 
stages of the EoR, when fluctuations in x H i are becoming influen¬ 
tial. However, even in this model the influence of x m is small, so 
the amplitude and position of the variance’s maximum should pro¬ 
vide a constraint on the X-ray production efficiency. The extent of 
the plateau that precedes it could provide insight into the relative 
timing between the onset of WF-coupling and X-ray heating. 

We can gain insight into the variance’s strong dependence on 
the correlation between Tf 1 and 5 by calculating the variance of 
A p,. We find that. 



and see that the variance is only sensitive to (T C mb/T b ) in the fi- 
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Figure 5. Redshift evolution of the brightness-temperature skewness (top) 
which exhibits a minimum associated with the onset of WF coupling, fol¬ 
lowed by a maximum that is driven by (T c mb /T s ) (bottom). Thin lines on 
the bottom plot correspond to the evolution of the ionized fraction for each 
model. Statistics are calculated from maps resolved to 4 Mpc, and a = 1.5 
for all models. 


nal term where its influence will be suppressed by a factor of 

{T CMB 5/T a ). 

In contrast, we find that the skew of A fi (for which we include 
the full equation in Appendix[A} is sensitive to both of these terms 
independently and in combination. The position and amplitude of 
the maximum in the skewness during this heating phase is mainly 
sensitive to (T C mb/T b ) as this factor dominates over (2cmb<5 /T b ). 
This is clear from Fig. [5] where we plot the skewness (top) and 
(Tcmb/Ts) (bottom) as functions of redshift. Initially the skew¬ 
ness becomes increasingly negative during the early stages of WF- 
coupling. There is a universal minimum to the skewness of our 
models at z ~ 31 driven by fluctuations in the WF-coupling (the 
details of which are unchanged between models) drawing the spin 
temperature towards the lower kinetic temperature (see the discus¬ 
sion surrounding Fig. |TJ. The skewness increases from this mini¬ 
mum, becoming positive and reaching a maximum as the average 
spin temperature (depicted in the top plot of Fig. [3} reaches its low¬ 
est point. At this point, the ft, parameter will be greatest and so 
fluctuations in the spin temperature dominate. 

As previously discussed, we see from the plot of (T C mb/T b ) 
in the bottom plot of Fig. [5]that the amplitude of the X-ray heating 
skewness maximum is inversely proportional to that of (T CM b/T b ). 
We find this to be due to contributions from negative ( T CM b/T b } 3 
terms becoming more dominant as the spin temperature decreases 
(see Appendix |A). 

Note that in the ( x < 10 56 models, the ionization field is be¬ 
coming influential as the skewness reaches its global maximum. If 
we plot the redshift evolution of {T (:MR x m /T B ) then we find a per¬ 
fect correlation between the peak in skewness and the minimum of 
this cross average. Even in such models, the high redshift maximum 
of the brightness-temperature skewness should provide constraints 


on the point at which the spin temperature is minimum, and thus 
the efficiency of X-ray production. 

|Shimabukuro et ah] < (2014| ) show the brightness-temperature 
variance and skewness for their fiducial model (£ x = 10 56 ). 
We mostly agree with their findings; however, their plot of the 
brightness-temperature variance only exhibits the X-ray heating 
peak (note their plot does not show the redshifts associated 
with reionization). The peak we associate with WF coupling and 
the plateau connecting it to the X-ray heating peak is totally 
absent. This may be because their boxes are small compared 
with ours. However, it is most likely that this difference is be¬ 
cause [Shimabukuro^eFar] ([20T4]( do not smooth their brightness- 
temperature maps prior to measuring one-point statistics, while we 
do(_] 


3.1.2 Hardness of the X-ray SED 

Fig. 0 (top) shows the redshift evolution of the brightness- 
temperature variance for different choices of spectral index, with 
= 10° 7 . The variance for the a = 3.0 (soft) model is more than 
double that of the a = 0.8 (hard) model. The softer the X-ray spec¬ 
trum the greater the anti-correlation between the density field and 
Tf (i.e. the spin temperature is smallest in underdense regions). 
This is evident in the bottom of Fig. [6] where we plot the redshift 
evolution of (T C mbS/T b ). This is to be expected as soft X-rays have 
a shorter mean free path than hard X-rays. 

The sensitivity of the variance amplitude to the spectral in¬ 
dex is degenerate with changes in amplitude produced by different 
X-ray efficiencies. This degeneracy maybe broken as the location 
and amplitude of the skewness’ X-ray heating peak is insensitive to 
variations of the spectral index (as seen in Fig. [7]in which we plot 
the skewness for different spectral indices, with £ x = 10 5 '). This 
is because, the redshift at which the spin temperature minimizes, 
and the difference between it and T CM b, is driven primarily by the 
efficiency of X-ray production. 

We expect the insensitivity of the skewness to the X-ray spec¬ 
tral hardness to be relatively model independent across the models 
we consider, as (T CM b 5/T s ) <C (T C mb/T s ) for all (see the bottom 
of Fig. [4] and Fig. [5}. However, should the X-ray production be so 
efficient that (T CM b/T b ) remains very small during this phase, then 
the skewness would be sensitive to (T C mbS/T s ), and therefore the 
X-ray spectral hardness. We conclude that if the efficiency can be 
constrained using the skewness, then the amplitude of the variance 
has potential for constraining the spectral index of the X-ray SED. 

|Pacucci et al.| < |2014] l find the peak amplitude of the large-scale 
(fc ~ 0.2Mpc~ ) power spectrum to be sensitive to the X-ray 
SED’s spectral index, but not the efficiency of X-ray production. 
We do not recover this behaviour by measuring the variance from 
maps smoothed on large scales. We find instead that, for smooth¬ 
ing scales of order 60 Mpc, sensitivity to the spectral hardness is 
lost whilst the amplitude remains sensitive to the efficiency of X- 
ray production. This suggests that the variance and power spectrum 
may be complementary in that the large-scale power spectrum can 
inform us on the spectral index and the variance smoothed on large 


18 There is a discretisation effect in 21CMFAST, associated with the gen¬ 
eration of the non-linear density field, that must be smoothed out in order to 
get a clean measure of the brightness-temperature statistics |Watkinson &| 
Pritchard 2014). This does not impact spin-temperature simulations, which 
are the focus of Shimabukuro et al~j^2i) 1 4^ 
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Figure 6. Evolution of variance with redshift (top) for various values of the 
spectral index of the X-ray SED; £ x = 10 57 for all. We find that suppres¬ 
sion of the (<5(Tcmb/7s)} amplitude caused by a harder X-ray spectrum, 
and seen in the bottom plot (for which the model key of the top plot ap¬ 
plies), reduces the variance. Thin lines on the bottom plot correspond to the 
evolution of the ionized fraction for each model. Statistics are calculated 
from maps resolved to 4 Mpc. 


Figure 7. Evolution of the brightness-temperature skewness with redshift 
(top) for various values of the spectral index of the assumed X-ray SED; 
£x = 10 57 for all. As (T C mb/T s — ) (shown in the bottom plot as a func¬ 
tion of redshift) is insensitive to the hardness of the X-ray spectrum, the 
skewness is not sensitive to the X-ray spectral index. Thin lines on the bot¬ 
tom plot correspond to the evolution of the ionized fraction for each model. 
Statistics are calculated from maps resolved to 4 Mpc. 


scales can provide constraints on the efficiency of X-ray produc¬ 
tion. 

3.2 Epoch of reionization 

When simulating reionization, it is often assumed that the spin 
temperature is totally saturated and therefore its fluctuations can 
be ignored. We see in Fig. [8] (in which we plot the brightness- 
temperature moments as a function of ionized fraction) that this 
may not be an appropriate assumption. Note that |Mesinger et ak| 
( [2013 1 discuss trends in the power spectrum's evolution at k = 
0.1 Mpc -1 similar to those seen in the variance. 


is more advanced when driven just by UV radiation. Such a shift is 
also seen in the minimum of the skewness and dimensional skew- 
ness^jassociated with ®i on ~ 0.2 (see the middle and bottom plots 
of Fig7[8] respectively). 

The late-time features of both skewness statistics at aiion ~ 
0.95, i.e. the rapid increase in the skewness as reionization ad¬ 
vances, and a turnover in the dimensional skewness, are far more 
robust. Although, for the highest efficiency we consider (£ x = 
10 58 ) the late-time turnover in the dimensional skewness doesn't 
occur, as X-ray ionizations complete reionization early relative to a 
UV-only model. In the middle plot of Fig. [3] we see that reioniza¬ 
tion completes at 2 ~ 9 in the £ x = 10 58 model; however, models 
which are mostly driven by UV ionizations don’t reach ii on ~ 0.95 
until a ~ 7.5. 


3.2.1 The impact of X-ray ionizations during reionization 

The variance (Fig.[8]top) for all our models is suppressed relative to 
that of a simulation that uses identical initial conditions but ignores 
spin temperature fluctuations (labelled here as ‘T s saturated’). This 
is due to partial ionization of neutral regions by X-rays. X-ray ion¬ 
izations are effective in both over and under-dense regions, reduc¬ 
ing the anti-correlation between the density and neutral-fraction 
fields. This is seen in Fig. [9] in which we plot (x m 8) as a function 
of ionized fraction, as the anti-correlation reduces with increasing 
X-ray efficiency. Partial ionizations also shift the mid-point maxi- 
murrrjin the variance to higher ionized fractions, as reionization 


3.2.2 The impact of heating on interpreting signatures of 
reionization 

Following the X-ray heating-dominated phase, discussed in Section 
m we see a rapid drop from the X-ray heating peak at low ionized 
fractions (in agreement with the findings of |Mesinger et al.|2013| l. 
If X-ray heating occurs relatively late as in the £ x = lO 88 model, 
the impact is dramatic as the drop from the heating peak occurs 


shows the PDFs during this phase. Unlike the £ x = 10 5 ' model 


during UV-driven reionization (occurring at Zion < 0.3)nFig. 


10 


19 The mid-point maximum refers to a maximum in the evolution of the 
variance during reionization. This occurs as the average ionized fraction 
of the Universe reaches 50% when the spin temperature is assumed to be 
saturated. 


20 the dimensional skew ness refers to th e skew normal ised with a 2 rather 
than cr 3 , this was found by |Watkinson & Pritchard|2014| to be a more natural 
ch oice during reionizatio n. 

21 |Mesinger et al. 1(2014) find that if X-ray heating is late enough, the heat¬ 
ing and reionization peaks can be merged into a single peak. 
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Average ionized fraction 

Figure 8. Brightness-temperature variance (top), skewness (middle) and di¬ 
mensional skewness (bottom) as a function of ionized fraction to highlight 
features relevant to reionization. X-ray ionizations are seen to reduce the 
amplitude of the variance (see the green dashed line with circles) by re¬ 
ducing the anti-correlation between density and neutral-fraction fields. If 
heating is very late (see the blue-dashed line) the turnover in the variance is 
no longer at the mid-point of reionization, which could lead to misinterpre¬ 
tation of the timing of the EoR. The skewness is qualitatively different in 
such a model, exhibiting a local maxima during the early stages of reioniza¬ 
tion (note that the nature of this feature will be extremely model dependent 
in this regime); such a feature could be used to constrain late X-ray heating 
models. Statistics are calculated from maps resolved to 4Mpc, and a = 1.5 
for all models. 

(where there is a clear distinction between a positive brightness- 
temperature distribution and a sharp spike at 5Tb = 0), the 
brightness-temperature distribution of the £ x = 10 5r ' PDF extends 
to negative temperatures. As a result, the contributions of neutral 
and ionized regions to the PDF are no longer distinct in brightness 
temperature. This reduces the variance and alters the skewness evo¬ 
lution, which exhibits a local maxima as the skewness tends to zero 
when 5Tb ~ 0 in neutral regions. 

Such signatures provide an opportunity to constrain the nature 
and timing of X-ray heating, fiowever they also complicate inter¬ 
pretation of the variance and skewness during reionization, impact¬ 
ing our ability to constrain reionization using these moments. For 



Ionized fraction 

Figure 9. Evolution of the (i h i<5) cross average as a function of ionized 
fraction; this confirms that it is the decrease in the negativity of (fern) that 
reduces the variance in the f x = 10 58 model. Statistics are calculated from 
maps resolved to 4Mpc, and a = 1.5 for all models. 
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Figure 10. Brightness-temperature PDF for redshifts (corresponding to 
Sion =[0.17,0.28,0.43]) associated with the local maximum feature seen 
in the skewness of the £ x = 10 55 model for: £ x = 10 57 in the top plot 
(where the zero peak corresponding to ionized pixels is separate from, and 
to the left of, the contribution from neutral regions); and £ x = 10 55 in 
the bottom plot (where the contribution from neutral regions is merged with 
that of the zero peak). Statistics are calculated from maps resolved to 4 Mpc, 
and a = 1.5 for all models. 

example, jPatil et al. |j2014} fit a function with a single peak to the 
variance of mock data, in order to constrain parameters of reioniza¬ 
tion. Such an approach would return misleading constraints, espe¬ 
cially if late X-ray heating occurred. 

|Ghara et al.| ( |2015^ note this fact and suggest to use either 
a three peak model (to model reionization, heating and coupling 
peaks) or a redshift cut-off. A redshift cut-off requires either prior 
knowledge on the timing of heating and reionization and/or throw¬ 
ing away information. The data itself could be used to provide a 
prior on where a redshift cut should be made (for example, model 
selection could be used to infer whether a three, two or one peak 
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model best describes the data in hand and where the transitions 
from one to the next occur). However, this would still be misleading 
in the £ x = 10 55 model, as the drop from the heating peak occurs 
over redshifts for which the ionized fraction is ~ 0.25 and the peak 
(usually associated with the mid point of reionization) is at ion¬ 
ized fractions of between 0.6 and 0.7. We therefore conclude that it 
would be prudent to use a parameter estimation approach that uses 
simulations to capture such subtleties. Unfortunately, this is par¬ 
ticularly challenging as simulations that include spin-temperature 
fluctuations are computationally expensive. Similar considerations 
would be necessary in constructing models of the skewness along 
the lines of jPatil et al.| ( [2014] >. 

These arguments are also relevant to MCMC parameter esti¬ 
mation using simulations that assume the spin temperature is satu¬ 
rated, such as that of Greig & Mesinger (2015) (who consider the 
power spectrum rather than the variance). It would be interesting to 
test the code they describe (21CMMC) against a mock dataset gen¬ 
erated from models similar to those we describe here to quantify the 
potential bias we would suffer from ignoring the spin temperature 
in performing parameter estimation. 


4 OBSERVATIONAL PROSPECTS 


To consider the effect of instrumental noise and foreground^ we 
make use of the publicly available code 21CMSENSE 22 (Pober et al. 


|2014|[20T3t . We refer the readers to the 21CMSENSE literature for 
details, but we will describe the main points for completeness. 

There are two main contributions to the error on the power 
spectrum: thermal noise and sample variance. At lower redshifts 
shot noise of the distribution of HI must also be considered, but 
this term is neglected in this analysis as it is found to be a sub¬ 
dominant effect, even after reionization (see |Pober et al.||2013) . 
|Pober et~af (2013 2014) calculate the noise for the fc-mode mea¬ 
sured by each individual baseline. As such, for a given redshift, the 
power spectrum may be calculated by application of an inverse- 
variance-weighted summation, for which the optimal estimator of 
the total noise error is. 


SAflk) = Y, 


[A* (fc) + A* (fc)]» 


( 12 ) 


The thermal noise contribution for a fc-mode labelled by i is given 
by A= XfYi[kf/(2Tv 2 )][Tli/(2ti)]Tf B . In this expression 
X 2 Y [h -1 Mpc] translates observed units into cosmological dis¬ 
tances; Cli [steradians] is the solid angle of the primary beam for 
a given baseline; TsysfK] is the system temperature, a combination 
of the sky and receiver temperatures, (i.e. T sya = T lec + T a k y ); and 
ti [hours] is the integration time for a given mode. The effect of 
the Earth’s rotation (relevant to a drift-scan observation modj^} 
is taken into account when calculating the noise on an individ¬ 
ual mode; i.e. different baselines may observe the same mode at 
different times which increases the integration time and therefore 
reduces thermal noise (for a similar reason, redundant baselines 
are useful). The sample variance error A 2 vi (fc) is equivalent to 
the 21-cm power spectrum for that mode at a given redshift, i.e. 
As v,i(k) — Ali,i(fc) = k 3 /( 2 n 2 )P 2 i(k) where P 2 i(k) is the 21- 
cm brightness-temperature power spectrum. 

As well as needing to beat down this error term, there is also 


22 https://github.com/jpober/21cmSense 

23 Note that SKA and LOEAR can also perform tracked scans. 


the issue of foregrounds which swamp the signal by several or¬ 
ders of magnitude. By considering the Fourier transform along the 
frequency axis of each mode independently (effectively measuring 
the delay in signal arrival time between the two interferometer ele¬ 
ments that make up a baseline), |Parso ns et al. ( |2012| ) find that the 
spectrally smooth nature of foregrounds mean that their contribu¬ 
tion will be confined to the region of delay space containing the 
maximum delays for a given baseline (confining them to be below 
an ‘horizon limit’). On the other hand, the 21-cm signal should ex¬ 
hibit unsmooth spectral characteristics so that some contribution 
from the cosmological signal will be observed with smaller de¬ 
lays (i.e. above the ‘horizon limit’). This motivates the definition 
of fcii min and k± :ln i n below which foregrounds will dominate. Be¬ 
cause of the frequency dependence of interferometer baselines, this 
‘horizon limit’ drifts to increasing values of fc|| with baseline length 
(i.e. with increasing k ±) producing a ‘wedge’ of foreground con¬ 
tamination in fc|| imin -fc^, m in parameter space. Mathematically the 
fc||,min ‘horizon limit’ may be described as i jParsons et al.|2012| , 

fc||,hor = , (13) 


where X and Y convert angle and frequency to co-moving distance 
respectively. 

There are two main approaches to dealing with the problem 
of foregrounds. One approach is to exploit the confinement of fore¬ 
grounds to the ‘horizon limits’ described above and essentially ig¬ 
nore modes that fall outside of EoR window (the region of fcii- 


k± space bounded by the ‘horizon limits’); see 

Datta et al. 2010 

Vedantham et al. 2012 Morales et al. 2012] f 

’hyagarajan et al. 

2013THazelton et al. 2013; Liu et al. 2014b When performing an 


inverse-variance-weighted (IVW) summation over fc-modes, this is 
equivalent to assigning infinite noise to modes that fall outside the 
EoR window. In parallel, there is a great deal of effort going into 
actively removing foregrounds from observations; these exploit the 
smooth spectral characteristics of foregrounds to identify and re¬ 
move their contribution (see Wang et al. 2006^ Liu & Tegmark 


201 l||Paciga et al.|201 l||Petrovic & Oh|201 1| Chapman et al.|2012[ 
Cho et al.|2012[|Shaw et al.|2014| >. 

Although the effectiveness of foreground removal has yet to 
be proved (for example, the impact of the frequency dependent na¬ 
ture of the instrument on the effectiveness of these removal tech¬ 
niques has yet to be established), we consider optimistically that 
it will be possible to remove foregrounds (described by ‘remove’ 
in the plots of IVW-brightness-temperature variance as a function 
of redshift in Fig. HD> and so reduce the wedge’s extent to the 
edges of the instrument’s field of view. In considering foreground 
avoidance (described as ‘avoid’ in the plots of Fig. |12| >, we as¬ 
sume that the spectral structure of the foregrounds only extend by 
Afc|| = 0.1/iMpc^ 1 beyond the wedge described by Equation 13 
(in line with the predictions of |Parsons et al.|2012} . For both fore¬ 
ground models we assume that baselines sampling the same k± can 
be combined coherently. 

We perform an IVW summation over the power spectra 
measured by an instrument (using errors from 21CMSENSE and 
the instrumental properties described in table to get con¬ 
straints on the brightness-temperature variancepj The 1 -<j er¬ 
ror on the IVW variance is estimated by Ei(l/^A 2 ) 2 ] -5 . 


25 In performing an inverse-variance-weighted summation over the power 
spectrum to calculate the IVW brightness-temperature variance we do not 
worry about normalisation factors as the power spectrum as a function of k 
is not bounded. As such, care must be taken in comparing the amplitude of 
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Figure 11. Inverse-variance weighted variance as measured when foregrounds are removed. LOFAR (Top-left); MWA (top-right); HERA (bottom-left); SKA 
(bottom-right). In the interest of brevity we do not explicitly consider PAPER’S performance, but PAPER-128 should be comparable in sensitivity to MWA 
for measuring the inverse-variance weighted variance. In this best case scenario all instruments are seen to be capable of constraining reionization; HERA 
and SKA would also tightly constrain WF-coupling and X-ray heating. 


Table 1. Summary of the instrumental properties we assume in calculat¬ 
ing our errors. We assume a bandwidth of 8 MHz and an integration time of 
1000 hours for all instruments. We use exact station locations for MWA and 
LOFAR from |Beardsley et ar]j2012) and jvan Haarlem et al.H2013) respec¬ 
tively (we also assume that each LOFAR HBA substation can be correlated 
independently to maximize redundant baselines). For HERA we assume 
331 antenna in a closely packed hexagon as per|Pober et al.|^2013f. Fol¬ 
lowing [PobeFeFar]J20T3J and |Greig & MesingerH2015) we model SKA 
as a closely packed hexagon (to maximize redundancy) of 218 antenna out 
to ~ 300 m with a further 215 stations randomly distributed in an annulus 
from to 300-600 m, 217 randomly placed stations in another annulus from 
600-1000 m. and 216 randomly placed in an annulus covering 1000-2000 m 
from the centre of the array. 


Parameter 

LOFAR 

MWA 

HERA 

SKA-1 

24 

Number of stations 

48 

128 

331 

866 


Element size [m] 

30.8 

6 

14 

35/v/2 

Collecting area [m 2 ] 

35,762 

896 

50,953 

416,595 

Receiver T [K] 

140 

440 

100 

40 



An inverse-variance-weighted sum over the dimensionless power 
spectrum (A 2 ) at different ki is only an unbiased estimator 
if the noise is Gaussian distributed with zero mean and the 
power is approximately flat (i.e. A 2 = A 2 for all i) so that 
(E I (A 2 + n ! )(l/5A 2 ) 2 /E 1 (l/5A 2 ) 2 } = A 2 . Of course this 

such an observation with simulations, i.e. the IVW variance must be simu¬ 
lated with equivalent modes to those probed by observations. 


is not strictly true as there is important evolution in the shape 
of the power spectrum with k. As such calculating of vw = 
Ei A 2 (l/<5A 2 ) 2 /Ei(l/5A 2 ) 2 means that of vw is sensitive to 
the details of the noise, which is in turn sensitive to the details of 
the instrument0The limited resolution of the instruments means 
that the power contribution from large k is totally suppressed; it 
is this that recovers the characteristics of the variance seen in our 
of vw statistic. Similarly, the IVW-variance is sensitive to the fore¬ 
ground model we assume; the presence of a foreground-corrupted 
wedge means that power from the associated k modes will be sup¬ 
pressed in calculating of vw - As is clear from the differing ampli¬ 
tudes between the plots of Fig. [TT] and those of Fig. [12] the level 
of foreground corruption can seriously impact the amplitude of the 
variance. We must therefore be very careful when interpreting this 
statistic quantitatively from observations. The power spectra mea¬ 
sured from fc-modes inside the EoR window should not suffer from 
this issue. The qualitative nature of the variance’s evolution is in¬ 
sensitive to the foreground corruption we consider here and could 
therefore be useful for constraining coupling, X-ray heating and 
reionization. 

We find that the first-generation instruments will only be able 
to constrain our models using the variance if foreground removal is 
possible. If so, then as is clear from the top row of plots in Fig. El 
both LOFAR and MWA will be able to constrain reionization and 
would also be sensitive to late X-ray heating. However, using fore¬ 
ground avoidance LOFAR could be sensitive to models in which 

26 Under the assumption that A 2 = A 2 for all i, then 
E ! A 2 (1/5A 2 ) 2 /E i (l/-5A 2 ) 2 = A 2 . 
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Figure 12. Inverse-variance weighted variance as measured from the 
foreground-free EoR window (i.e. by assigning infinite noise to power from 
foreground corrupted modes) by LOFAR (top), as well as by HERA and 
SKA (bottom). In this worst case scenario, LOFAR will be limited to plac¬ 
ing upper bounds on the EoR, while HERA and SKA could still tightly 
constrain the EoR. HERA and SKA will still be able to constrain X-ray 
heating, but not WF-coupling. 


reionization ends later than our models assume, but will more likely 
be limited to setting upper limits (see the top plot of Fig.[T2]l. Note 
that because of the maximal redundancy of its baselines, the next 
phase of PAPER (consisting of 128-elements, see |Ali et al.fl20l5) 
is only marginally less sensitive than MWA (see the appendix of 
|Pober et al.|2014} despite having less than half the collecting area. 

We again emphasise that our fiducial reionization model is op¬ 
timistic, and so first-generation instruments may struggle more than 
is suggested by our analysis. Furthermore, due to the presence of 
sinks in the IGM the variance may be up to a factor of two smaller 
than the models of this paper predict. If extreme levels of rem¬ 
nant HI in galaxies are present then the variance will be reduced 
even further (Watkinson et al.||2015| . Such reduction of the vari¬ 
ance could make it difficult for first-generation instruments to do 
more than place upper limits, even if foregrounds can be removed. 
However, the IVW-variance will inevitably have smaller errors than 
the power spectrum at a given k and therefore, first-generation in¬ 
struments would do well to exploit it in their quest to make a first 
detection. 

Under the same assumptions, next-generation instruments 
such as SKA and HERA will be able to tightly constrain the vari¬ 
ance for the coupling, heating and reionization epochs (as seen in 
the bottom plots of Fig. HQ- Note that the IVW variance exhibits 
three distinct peaks corresponding to WF coupling, X-ray heating, 
and reionization; it does not exhibit the plateau between WF cou¬ 
pling and X-ray heating as seen in the standard variance (i.e. that 
measured from a clean simulated box). Even if foreground removal 


proves intractable then the foreground-avoidance technique will re¬ 
turn strong constraints on the heating and reionization epochs (see 
the bottom plot of Fig. m Sinks in the IGM will not stop these 
next-generation instruments from returning strong constraints on 
the EoR using the moments. However, extreme levels of residual 
HI can qualitatively alter the evolution of the moments from that 
described in this paper (Watkinson et al. 20151. 

We use the approach detailed in Watkinson & Pritchard| l [2014| > 
to approximate instrumental errors on the skewness, this approach 
assumes that foregrounds can be perfectly removed and approxi¬ 
mates instrumentals by smoothing and re-sampling pixels to match 
the resolution of the telescope. We plot the dimensional ske wness 
as a function of both ionized fraction and redshift in Fig. 


13 


W 


These errors should be viewed as optimistic estimates and will 
likely be quite a bit larger. As an illustration, if we compare the er¬ 
rors on the variance as calculated by |Watkinson & Pritchard||2014| > 
with those predicted for the IVW variance, we find its S/N is a fac¬ 
tor of order 3 worse if foregrounds can be removed; if foreground 
avoidance is necessary then S/N can be 20 - 50 times worse. 

We see that it will be possible to use the skewness to constrain 
models of late X-ray heating, possibly with LOFAR but certainly 
with the next-generation instruments. Therefore this presents an ex¬ 
cellent opportunity for these telescopes to constrain a fundamental 
property of the Universe’s evolution, namely the relative timing of 
WF coupling. X-ray heating, and reionization. 


5 DISCUSSION 

There are several approximations made in 21CMFAST that may 
have important repercussions; in particular the code assumes aver¬ 
age properties of the IGM in calculating the X-ray mean free path. 
This is most important during the later stages of reionization when 
large ionized regions will result in fluctuations of the X-ray mean 
free path between sight-lines in the box. The code also assumes ei¬ 
ther population II or III stellar spectra in its calculations of WF cou¬ 
pling, and does not account for the possibility of mixed populations, 
feedback effects or shot noise. The nature of these first stars (and 
of the remnants they leave when they die) is very uncertain. Recent 
simulations indicate that these first stars will be 10 — 100 Af© and 
are expected to form in small clusters (e.g. |Hirano et al.|2014||Greif| 
|et al.|20fT| >. Formation of such population III stars rely primarily 
on cooling via molecular hydrogen, however they produce large 
amounts of Lyman-Werner radiation (which disassociate molecular 
hydrogen) and so are likely to stunt further formation of population 
III stars (e.g. |Wise & Abel|2007||Q’Shea & Norm an 2008). Such 
large stars are also short lived, so it is not unreasonable (as is done 
in this work) to assume that population II stars will be the dominant 
driver of the processes discussed here. However, it is possible that 
these results are inaccurate during the very early phases of coupling 
when the very first stars form. 

Whilst simulations such as 21CMFAST have been tested 
against numerical simulations during the EoR assuming that the 
spin temperature is saturated (see for example, [Zahn et al.[201 l| and 
|Majumdar et al.|2014| >, there has not been equivalent tests of these 
when spin temperature fluctuations are included. This is mainly 
because numerical simulations with the necessary scale and res¬ 
olution do not yet exist. The only numerical simulations that per¬ 
form radiative-transfer in all of the relevant frequency bands are 

2 ' Prior to measuring the skewness for this figure, we smoothed brightness- 
temperature boxes to a radius of lOMpc to suppress noise corruption. 
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Average ionized fraction 



Figure 13. Dimensional skewness as a function of ionized fraction (top) 
and redshift (bottom) from brightness-temperature maps smoothed to a ra¬ 
dius of lOMpc (top). Errors (shaded from dark to light for SKA, HERA 
and LOFAR respectively) correspond to approximate instrumental errors 
only and are therefore optimistic estimates. Whilst the first-generation in¬ 
struments will struggle to extract useful constraints from this statistic, the 
skewness measurements from HERA and SKA will return vital constraints 
on the EoR. When we examine this statistic as a function of redshift (bot¬ 
tom), we find that SKA will also be able to tightly constrain X-ray heating 
with the skewness, but that HERA will not be sensitive to the skewness 
beyond the EoR. 


those of |Baek et al.| l |2010| ). These simulation do not resolve haloes 
with mass below 10 1u Mq, therefore they do not resolve atomically 
cooling haloes. As such, all astrophysical processes are driven by 
more massive, and therefore more rare and biased, haloes than is to 
be expected in reality. It is therefore not possible to draw direct 
comparison between [Baek et aI7| < |2010| ) and 21CMFAST . How¬ 
ever, |Mesinger et al.| |20 1 3[ > note that the qualitative evolution of 
the power spectrum at k ~ 0.1 of 21CMFAST (when including 
spin temperature fluctuations) is in agreement with the numerical 
simulations of jBaek et al.| ( j2010) . We also find that the skewness of 
our late X-ray heating model (£ x = 10 55 ) qualitatively agrees with 
their S6 model, which is encouraging. 

There are other processes that must be considered in paral¬ 
lel to spin temperature fluctuations. For example, and as already 
discussed, the presence of sinks could drastically reduce the vari¬ 
ance. This reduction is due to residual signal in ionized regions and 
sub-pixel ionized regions. X-ray ionizations will occur in a more 
homogeneous fashion than UV ionizations and so will be responsi¬ 
ble for partially ionizing regions outside of UV carved ionized re¬ 
gions. It therefore seems likely that the reduction of variance caused 
by X-ray ionizations will be in addition to that caused by sinks, 
i.e. they will further reduce the contrast between over and under- 
dense regions. The simulation of |Sobacchi & Mesinger|j2014) also 
incorporate UVB feedback which suppresses star formation, such 


feedback will clearly impact on both the Lyman-a and X-ray pro¬ 
duction. However, given the large amplitudes seen in the one-point 
statistics during the heating epoch, and that we have studied four 
orders of magnitude in the X-ray efficiency, it is unlikely that UVB 
feedback will have a dramatic effect beyond that seen here. 

These examples (and the lack of numerical simulations with 
which to test 21CMFAST ) serve to illustrate the challenge we face 
in simulating the epochs of the first dawn and reionization. The 
results of this work should therefore not be considered conclusive 
and it is essential that we do more to understand how the statistics 
of the 21-cm moments are impacted by different physical processes 
(and their interplay). 


6 CONCLUSION 

In this paper, we have considered the sensitivity of one-point statis¬ 
tics of the 21-cm brightness temperature to fluctuations in WF cou¬ 
pling and X-ray heating, concentrating on the skewness and the 
variance. We use semi-numerical simulations to vary the efficiency 
at which X-rays are produced (to cover four orders of magnitude) 
and the spectral index of the X-ray SED (to encompass the range of 
observational constraints we have at low redshifts). From this study 
we establish that: 

(i) the location and amplitude of the global maxima in the red¬ 
shift evolution of both the skewness and variance are sensitive to 
the X-ray production efficiency. The amplitude of this maximum in 
the variance is also sensitive to the hardness of the X-ray SED. This 
degeneracy may be broken, as the skewness is only sensitive to the 
X-ray production efficiency; 

(ii) late X-ray heating causes the drop from the X-ray heating 
peak to occur at an ionized fraction of about a quarter rather than in 
the very early stages of reionization. In such a model, the turnover 
in the variance, usually associated with the mid-point of reioniza¬ 
tion, is shifted to higher ionized fractions. The evolution of the 
skewness is qualitatively different if X-ray heating occurs late, this 
provides a clean way to constrain such a model. The amplitude 
of the variance is greatly reduced in these models, which would 
make it more challenging for the first-generation instruments (such 
as LOFAR, MWA and PAPER) to make a detection of reionization 
using the variance; 

(iii) the high-redshift heating peak must be allowed for in mod¬ 
els used for parameter estimation from one-point statistics. If not 
our inferences may be very misleading. This is equally true for per¬ 
forming parameter estimation from the power spectrum; 

(iv) X-ray ionizations reduce the amplitude of the variance. In 
most models we consider they reduce the variance by ~ 10% 
during the mid-phase of reionization; in the most X-ray efficient 
model, we find this reduction to be ~ 25%. 

We consider (for the first time to the authors’ knowledge) the 
variance as measured using foreground avoidance techniques. From 
this we find that the next-generation instruments such as HERA and 
SKA will return strong constraints on both reionization and X-ray 
heating, even if we are unable to remove foregrounds. 

The findings of this paper will help us to correctly interpret 
future observations of the 21-cm brightness temperature; in par¬ 
ticular they have important consequences for improving parameter 
estimation during reionization. 
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APPENDIX A: ANALYTICAL EXPRESSION FOR THE 
SKEW OF A fj, 
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